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We develop the formalism for the evaluation of density-density correlators in lattice QCD that 
includes techniques for the computation of the all-to-all propagators involved. A novel technique 
in this context is the implementation of the one-end trick in the meson sector. Density-density 
correlators provide a gauge invariant definition for the hadron wave function and yield information 
on hadron deformation. We evaluate density-density correlators using two degenerate flavors of 
dynamical Wilson fermions for the pion, the rho-meson, the nucleon and the A. Using the one-end 
trick we obtain results that clearly show deformation of the rho-meson. 



I. INTRODUCTION 

Deformation in nuclei [TJ |2] and atoms [3J 3] is an im- 
portant phenomenon that has been extensively studied. 
In this work we address the question of whether defor- 
mation also arises in low-lying hadrons using the funda- 
mental theory of the strong interactions, Quantum Chro- 
modynamics defined on the lattice. In order to be able 
to answer this question we develop techniques for the ex- 
act evaluation of four-point correlators. These methods 
are also needed in a range of other applications in lattice 
QCD. 

In this work we study the shape of the pion, the rho- 
meson, the nucleon (N) and the A. The pion being a 
spin-0 particle is expected to have no deformation and it 
therefore provides a check for our methodology. For par- 
ticles with spin larger than 1/2, the one-body quadrupole 
operator provides a convenient characterization of defor- 
mation. The spin 1/2 nucleon cannot have a spectro- 
scopic quadrupole moment but can still have an intrin- 
sic deformation. The experiment of choice to reveal the 
presence of deformation in the nucleon and its excited 
state the A is measuring the N to A transition ampli- 
tude. Significant effort has been devoted to photo- and 
electro-production experiments on the nucleon at major 
experimental facilities [3 O [3 |5]. These experiments 
measure to high accuracy the ratios of the electric (E2) 
and Coulomb (C2) quadrupole amplitudes to the mag- 
netic dipole (Ml) amplitude. If both the nucleon and 
the A are spherical, then E2 and C2 are expected to be 
zero. There is mounting experimental evidence over a 
range of momentum transfers that E2 and C2 are non- 
zero [9, 10 . These ratios have been recently shown to be 
non-zero in lattice QCD |llj pointing to deformation in 
the nucleon or/ and A. 

A different approach that sheds light on deformation 
is to use density-density correlators to directly probe the 
hadron wave function [T2 US] . Density-density correla- 
tors [lllllSlIlSlITllllHlIiniEOlETlEa provide 
a gauge invariant definition of the hadron wave func- 
tion. In a previous study |16j the density-density cor- 
relators were evaluated approximately. This was due to 
the fact that the all-to-all propagators needed for their 
exact evaluation were not calculated. Furthermore they 



were computed for pion masses larger than 600 MeV and 
on lattices with a spatial volume of about 1.5 fm. 

In this work we provide an exact evaluation of the 
four-point functions involved in the computation of the 
density-density correlators. The all-to-all propagators 
needed for the exact evaluation are calculated using 
stochastic techniques combined with dilution. In addi- 
tion, we apply in the meson-sector for the first time in 
this context, the so-called one-end trick originally de- 
vised to evaluate the pion zero momentum two-point 
function |25j . In the two-point function, the one-end 
trick amounts to a clever summation of the spatial co- 
ordinates not only of the sink as routinely done but also 
of the source and therefore all-to-all propagators are in- 
volved. Implementation of this trick in the evaluation of 
the meson density-density correlators leads to a signifi- 
cant reduction of the statistical errors [24]. This trick, 
in its present formulation, can only be applied to meson 
density-density correlators. In baryons, the density in- 
sertions are on only two of the three quarks which gives 
rise to an odd number of quark propagators that cannot 
be grouped in pairs for the summation to work. 

An alternative method applicable to both mesons and 
baryons is to combine stochastic evaluation of one all- 
to-all propagator with a sequential inversion to sum over 
the other spatial coordinate. This method, apart from 
the requirement of fixing the final hadronic state, needing 
new sequential inversions for each of the nucleon and A 
states, has been shown to yield results with similar errors 
as using two sets of stochastic inversions [22] We therefore 
do not consider it here. 

Further improvements as compared to the previous 
study of density-density correlators [TB] is that we use 
a spatial lattice of 2A^ as compared to 16"^ used previ- 
ously and dynamical Wilson fermions corresponding to 
smaller pion masses, the lowest being 380 MeV. 

This paper is organized as follows: In Section II we de- 
fine the density-density correlators, in Section III we ex- 
plain the stochastic techniques used for the evaluation of 
the all-to-all propagators, in Section IV we give the inter- 
polating fields and parameters of the simulations and in 
Section V we describe our results on the density-density 
correlators for the pion, the rho-meson, the nucleon and 
the A and show how to correct for finite spatial volume 
effects. Finally in Section VI we summarize and give our 
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conclusions. 



II. DENSITY - DENSITY CORRELATORS 

Throughout this work we consider the equal-time 
density-density correlators defined by: 



d^xi [ d^x{n\jH{x,t)j^{x2 + xi,h)j^{xi,h)jl{xo,to)\n) (1) 



where Jq is the normal ordered density operator : qjoq : and Jh is an interpolating field with the quantum numbers 
of the lowest lying hadron H. The two integrals in Eq. (jlj ensure that the state is projected to zero momentum; one 
integral sets the momentum of the sink equal to that of the source while the other sets both to zero. This can be 
shown explicitly by inserting three complete sets of states in Eq. ([ij : 

-£;„^.(0)(t-ti) ip-d-a g-B„,(0)(ti-to) 

CH{S2,ti)= (^\JH\nf,6) (n^,Ob-|n,p)— — (n,p1jg|n,,0) (n.»0|4|»). (2) 

p,n,n^,nf 2A„^ (0) ^^n(P) 2A„. (0) 

In the large ti — to and t ~ ti limit we have: 

— »-oo;(fi— io)— >oo CH{x2,tl) 

p-mff(*-to) pip-S2 

= Y^mjum^ ^^2^ {m\n,P) {n.mH). (s) 



If we divide by the zero momentum hadron two-point 
function G/f (0, t~to) then the exponential dependence on 
t — to and overlaps cancel and we obtain the expectation 
value of the two density insertions, |jg (a;2)jo 1-^)- 
the non-relativistic limit, this expectation value gives the 
charge distribution of the hadron. It can be written in 
terms of the non-relativistic form factors [13] 

pip-X2 

(^bo(^2)jo'l^> ^T.^HniP) FnHi-P) (4) 

p,n 

where 

F]^j„{p) ^ {H\j^\n,p) . (5) 

The connected diagrams of the density-density corre- 
lators for mesons and baryons are shown in Fig. [l] We 
note here that the diagram depicted in Fig.[l]for baryons 
yields a correlator that depends only on one relative dis- 
tance instead of two. To obtain, in the non-relativistic 
limit, the charge distribution that depends on the two 
relative distances one must calculate the three-density 
correlator. This requires the evaluation of two types of 
five-point functions shown in Fig. |2] In Ref. [161 the 
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FIG. 1: Equal-time density-density correlators for mesons 
(upper diagram) and for baryons (lower diagram). 



three-density correlator or five-point function was eval- 
uated approximately for one of the diagrams shown in 
Fig. [2] for which each quark line has only one density in- 
sertion. It was shown that integrating over one relative 
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distance one obtains results that are consistent with the 
corresponding two-density correlator. For the work pre- 
sented here we therefore only consider correlators with 
two-density insertions, which give the distribution of one 
quark relative to the other irrespective of the position of 
the third. In other words, in the non-relativistic limit, it 
corresponds to the one-body charge distribution. 




FIG. 2: The three density correlator for baryons. 



What makes four-point functions harder to evaluate 
than three-point functions is the fact that we need to 
compute all-to-all propagators. Sequential inversions 
used in the evaluation of three-point functions can not 
be used here. The reason is that we are interested in 
obtaining the dependence in terms of a relative distance 
and therefore the spatial positions where the density op- 
erators are inserted involve the relative distance and can 
not be summed independently. Therefore the bulk of this 
work deals with the evaluation of the all-to-all propaga- 
tors to sufficient accuracy. 



III. STOCHASTIC TECHNIQUES 

The technically challenging aspect of the calculation 
of the density-density correlators is the fact that the 
summation over sink and insertion coordinates requires 
knowledge of all-to-all propagators. A previous study 
has been carried out in the quenched approximation 
and using two dynamical degenerate Wilson fermions in 
which no summation was performed over the sink coor- 
dinates . This eliminated the need of calculating all- 
to-all propagators at the cost of not explicitly projecting 
to zero momentum states, which instead were only ob- 
tained via the large Euclidean time suppression of higher 
momenta. In this work we use stochastic techniques to 
estimate the all-to-all propagators f26', '27' enabling us to 
sum over the sink coordinate and thus explicitly project 
to zero momentum initial and final states. 

In order to evaluate the all-to-all propagator one be- 
gins by defining an ensemble of Nr noise vectors ^^(x, t)r 

obeying to order (y^) 

{^;,{x,t))r - and 

where /z and a are spinor and color indices respectively 
and r enumerates the vector in the stochastic ensem- 
ble. In particular, we use Z(2) noise where ^^(af, i) G 
— 1,— i} with equal probability. By solving the 
Dirac equation with each of these Nr noise vectors as 



the source, one obtains an ensemble of solution vectors: 
y 

where (f> is a, solution vector and G is the inverse of the 
Dirac operator. If we now take the average over the prod- 
uct between solution and noise vectors over the stochastic 
ensemble, we obtain an estimate of the all-to-all propa- 
gator: 

z 

Z 

= G;lix;y). (8) 

A well known technique used to suppress stochastic 
noise is dilution [28]. Within this technique, one dis- 
tributes the elements of a noise vector over certain color, 
spin and volume components of multiple noise vectors 
setting the remaining components to zero. An example 
is spin dilution where the first noise vector has non zero 
entries only on the first spin component, the second vec- 
tor only on the second spin component and so on. In 
this example, in order for the conditions in Eq. ^ to 
be satisfied, the total number of noise vectors in the 
ensemble is restricted to multiplets of four. In Fig. |3]we 
show a schematic representation of n-fold dilution. 
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FIG. 3: A schematic representation of n-fold dilution. Zi 
denotes a random complex number. 

The more one dilutes, the closer an estimate one gains 
of the all-to-all propagator. This can be understood if one 
considers the extreme case where a noise vector is diluted 
over all color, spin and volume components. In this case 
one would have inverted for each color, spin and volume 
index thus obtaining the exact all-to-all propagator. 

The straight forward way to carry out the computation 
of the density-density correlator is to expand Eq. (jlj on 
the quark level and replace each all-to-all propagator with 
the stochastic average over the product between solution 
and noise vectors: G'^j,{x;y) = (0^(a;)^J''(y)),-. Through- 
out this paper we will refer to this as the direct method. 
As demonstrated in Section IV, a reasonable estimate of 
the all-to-all propagators can be computed through the 
direct method if a large enough number of stochastic in- 
versions is carried out. 

Significant improvement to the results obtained using 
the direct method is achieved by applying the so called 
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one-end trick. The one-end trick was originally devised 
to compute pion two-point functions [25j . In its original 
form it is based on the realization that by appropriately 
combining solution vectors one can derive the pion two- 
point function summed over both ends (source and sink) . 
To be specific, let us consider the pion two-point function 
which, at the propagator level, is just the trace of the 
absolute square of the quark propagator: 



^(^(f , i)|7r(afo, to)) = ^ [\G{x, t; xo, to) H . (9) 



Let us consider the stochastic average over the product 
between two solution vectors given by: 



t; to)(j){x, t; to))r, 



(10) 



where the appearing in the argument of the solution 
vector is to indicate that the noise vectors are localized 
on this time slice, i.e: 



C{x,t)r^C{x)rS{t-to), 



(11) 



and hence 



^lix, t; to)r = J2 G*;^' y' to)Ciy)r- (12) 



By substituting for (/)° in Eq. (10) we obtain 



J2{'l^7is,t;to)K(^,t;to))r ^ 

X 

^ G;f{x; x'o)Gll{x; x'^)5,,5.J{x'o ~ 4') 
Y,Tr\\Gf^{x;x'o)\ 



(13) 



where Xg = {to^x'o) and x'o — {to^x'^). This is the pion 
two-point function given in Eq. ^ summed over all spa- 
tial source and sink coordinates. This double summation 
increases statistics by spatial volume as compared to the 
standard way where one computes two-point functions 
using a point-to-all propagator. The increase by spatial 
volume in statistics far outweighs the stochastic noise in- 
troduced by the stochastic inversion. 

The pion two-point function is the simplest implemen- 
tation of the one-end trick since the 7-structure of the 
interpolating fields combined with the backward propa- 
gator of the antiquark yield a simple trace over a prod- 
uct of two forward quark propagators. To apply the 
trick on an arbitrary meson two-point function with in- 
terpolating operators of the form qfTqf, where f =/= f 
label two flavors of quarks, not necessarily degenerate 
and r an arbitrary combination of gamma matrices, one 



must use spin dilution. More explicitly, the noise vec- 
tors should be of the form £,^{x)(^r,a) = ^°'{x)rSfj.a- The 
r index counts sets of noise vectors, each set contain- 
ing four noise vectors carrying an index a. We note 
here that this form of dilution is different than that de- 
scribed in the previous section. Here the Z(2) random 
numbers involved in the spin dilution are the same for 
each spin component entry. It can be easily confirmed 
that this choice satisfies the conditions in Eqs. (|6|; the 
sum over the stochastic ensemble now becomes a double 
sum (over r and a) and {^"'{x)^^"- {x'))r — 5{x — x')5aa'- 
Within this notation the solution vectors are denoted as 
'l^tii.x)(r,a) = Y.xa GfA^] xq)^^ {xo)r5ai, ■ Now One Can ap- 
propriately combine the solution vectors to incorporate 
the r matrices involved and obtain the meson two-point 
function summed over both ends: 

^ 0^(:?, t; to)(rM)^'va'l>T{x, t; to)(r,a)T'i^f, = 

x,r 

J2 Gf,{x- x'o)V'^,Glf {x- x'^)f',^5{x'o - x'^)5,,. = 
Y^Tt[G{x-xo)VG{xo;x)V] (14) 

where V = and f = 7or^7o. Thus the one-end trick 
can be generalized to an arbitrary meson interpolating 
field. We would like to note here that the automatic 
summation over the source using the same set of solution 
vectors selects a given momentum. Therefore the one-end 
trick by construction computes only two-point functions 
at a specific momentum. In the examples given above 
this momentum was set to zero. To compute meson two- 
point functions at various momenta, one must invert for a 
new set of solution vectors having previously transformed 
the noise vectors with an appropriate phase. In other 
words, one needs a new set of stochastic inversions for 
each momentum vector. 

The crucial point that makes the one-end trick applica- 
ble to the evaluation of density-density correlators is the 
fact that the initial and final states have zero momen- 
tum. To show how to implement the one-end trick we 
consider the density-density correlator for an arbitrary 
meson with an interpolating operator of the form qfTqf, 
where / 7^ /': 

C{x2) = ^Tr [757oG'(xi; xo)f 'G'''(x2+i; Xo) x 

757oG(x2+i;x)r'Gt(xi;x)] (15) 

where X2+1 = {ti,X2 + fi), xq = (^0,^0), xi = (ii,xi), 
X = [t, x) and V = T^^. Let us define: 

SfA^'^x;y\to) = ^0",(x;to)(r,CT)raKC^(y;io)(r,K) 

r 

(16) 

where x — {tx,x) and y — (ty,y) and the to appear- 
ing in the argument of 5"°^ is to indicate that the noise 
vectors are localized on time-slice in. Summation over all 
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repeated indices is implied. Assuming that the noise vec- 
tors are spin diluted in the manner described previously, 
we obtain 



5^^(r;x;y;io) = 

X! C!''^li^^tx;xo,to)T„^G*J'^ {y,ty;yo,to)6a'b'S{xo - yo) 
So-.yo 

= ^G{x,t^;xo,ta)T G\y,ty; xo,to)\'^^^ . (17) 



Thus in terms of the propagator defined in Eq. (16), the 
expression 

^ Tr [757oS'(f '; xi;x2+i;ta)j5'^oS{T'; X2+1; xi;t)] 



(18) 
15 ) with 



yields the density-density correlator of Eq. 
an additional summation over the source coordinate xq. 
This is the generalization of the one-end trick to meson 
four-point correlators. It is apparent that one needs two 
sets of stochastic inversions: one with the noise vectors 
localized on the source time-slice to and one with the 
noise vectors localized on the sink time-slice t. 



IV. INTERPOLATING FIELDS AND LATTICE 
PARAMETERS 

For the pion and the p-meson we compute the density- 
density correlators using both the one-end trick and the 
direct method. For the nucleon and the A it is not as 
straight forward to apply the one-end trick. The quark 
line propagating without a density insertion complicates 
the generalization of the trick to baryons since the propa- 
gators to be replaced by noise vectors are odd in number 
and therefore unlike for mesons the noise vectors cannot 
be grouped in pairs to yield (5-functions after summa- 
tion. Thus in this work for the nucleon and A density- 
density correlators we only present results using the di- 
rect method. 

One of our main goals is to detect a possible asymme- 
try in the charge distributions of these particles. For this 
purpose we select interpolating operators so that they 
project to physical spin states. For the mesons we use in- 
terpolating operators of the form: J*^ = uTd with F = 75 



for the case of the pion and F = { 



7l-»72 



71+172} for 



the -1-1, and —1 polarizations of the vector meson re- 
spectively, where we have taken the z axis along the spin 
axis. For the nucleon we use = e'''"^u°(u''^ 6*750?'^) 
where C — 7072- For the case of the A we opt to probe 
the spin ±| components. Thus we use the interpolating 
operators: 

J\ = 



1 ,a6c 



uJ(2M''^CF+d=) + d'j'(u''^CF+M=)] 



J^3 - 



[u|(2u''^CF_d^) + (u''^CF_u'=)] 



(19) 



where r± = (71 =F 172) /2. 

Given the large number of inversions needed to com- 
pute the density-density correlators and the available 
computer resources, using dynamical Wilson fermions 
that are fast to invert is the only option at our dis- 
posal. We use two dynamical degenerate flavors of Wil- 
son fermions at three pion masses. The exact parameters 
of the ensembles used are listed in Table HI 



TABLE I: The first column gives the number of configurations 
analyzed, the second the value of the hopping parameter, the 
third the pion mass in GeV, the fourth the ratio of the pion 
mass to the p mass, the fifth the nucleon mass in GeV and 
the last column the size of the lattice. The first two sets 
of configurations are from Ref. [29 while the third is from 
Ref. 130j. The lattice spacing is determined from the nucleon 
mass at the chiral limit. 

P = 5.6, a~' = 2.56(10) GeV 
iVconf K m^ (GeV) m^/mp Mat (GeV) L^xT 

185 0.1575 0.691(8) 0.701(9) 1.485(18) 24^x40 
150 0.1580 0.509(8) 0.566(12) 1.280(26) 24^x40 
200 0.15825 0.384(8) 0.453(27) 1.083(18) 24^x32 
= 0.1585 0.938(33) 



To suppress excited state contributions we use Gaus- 
sian or Wuppertal smeared sources [31]. In addition we 
apply hypercubic (HYP) smearing [37 on the gauge links 
that enter the smearing function that builds the Gaussian 
smearing function. The parameters that enter the Gaus- 
sian smearing function are taken from Ref. |33j where 
they were determined by optimizing ground state domi- 
nance for the nucleon. In fact, in Ref. |33j it was demon- 
strated that one can damp excited state contributions to 
the nucleon two-point function as early as 0.3 fm from the 
source time slice. The parameters for the HYP smearing 
are taken from Ref. ^2^. 

For the computation of the correlators we take the 
time-slice of the density insertions to be at mid-point 
of the time separation between sink and source. For the 
direct method we take the time separation between the 
sink and the source to he t — to — 10a or 0.77 fm. This 
is the minimum time separation that is needed for the 
suppression of excited states. For the one-end trick the 
separation between sink and source is set to i — to = 14a. 
The reason for taking a larger time separation when using 
the one-end trick lies in the accuracy of the results that 
allows for a larger time separation with a good signal. 
This allows us to check that indeed excited state contri- 
butions are sufficiently suppressed by comparing results 
at the two sink-source time separations. 

We first give the details of the computation in the case 
of the direct method. We require two sets of stochastic 
propagators per configuration, one with the noise vec- 
tors localized on the insertion time-slice and one with 
the noise vectors localized on the sink. We also compute 
a point-to-all propagator from the source time-slice to all 
lattice sites. The noise vectors are diluted in color, spin 
and even-odd spatial sites. Dilution in time is automatic 
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here since we invert with the noise vectors localized on 
a single time-slice. Thus each noise vector is diluted to 
twenty-four independent noise vectors requiring twenty- 
four times more inversions. The number of noise vectors 
used is determined through a tuning process. For this 
tuning the A-baryon correlator at the lightest pion mass 
is considered. By comparing the decrease of the relative 
statistical error when increasing on one hand statistics 
and on the other hand the number of noise vectors used, 
we determine the optimum number of stochastic vectors. 
For this tuning we use 50 configurations and compute the 
A-baryon correlator for three, six and nine such 24- fold 
diluted noise vectors. For Nr=3, 6 and 9 we find a rel- 
ative statistical error of 50%, 20% and 16% respectively. 
The fact that by doubling the number of noise vectors 
from 3 to 6 the statistical error decreases by more than 
one half is an indication that Nr=3 is too small yielding 
large stochastic noise. On the other hand, increasing the 
number of noise vectors from 6 to 9 the relative error 
decreases by a/6/9, which is what is expected from scal- 
ing. This indicates that at this point increasing Nr or the 
number of configurations is equivalent. We thus fix the 
number of noise vectors to six. Since we carry out two 
sets of stochastic inversions, one at the sink and one at 
the insertion time-slice, and since we use color, spin and 
even-odd dilution we need 288 stochastic inversions per 
configuration. This amounts to a total of 300 inversions 
per density-density correlator if we additionally consider 
the point-to-all propagating from the origin. To increase 
statistics for the two ensembles corresponding to the two 
hghtest pion masses needed for the baryons, we calculate 
density-density correlators using the first and second half 
time interval of each configuration. Furthermore, for the 
lightest pion mass we improve statistics by using N^—Q 
noise vectors for the correlators. Thus for k — 0.1580 
we carry out 600 inversions per configuration while for 
K = 0.15825 888 inversions per configuration. 

For the case of the mesons we have additionally com- 
puted the charge distributions using the one-end trick. 
Therefore for the computation of the meson density- 
density correlators additional inversions are carried out 
since the dilution method is specific to the one-end trick. 
Like for the direct method, two sets of inversions are 
needed to extract the density-density correlator using the 
one-end trick: one set with the noise source set at the 
source time to and one set with the noise source set at 
the sink time t. We use eight spin-diluted noise vectors 
amounting to 32 inversions at the source and 32 at the 
sink or a total of 64 per configuration. 



RESULTS 



Comparison between the direct method and the 
one-end trick 



For the meson density-density correlators we can com- 
pare results obtained using the direct method with those 



using the one-end trick. Given that the time separation 
between sink and source is larger in the latter case this 
also provides a check of ground state dominance. 
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FIG. 4: The pion density-density correlator using the one- 
end trick (upper graph) and using the direct method (lower 
graph). The mean value of C7r(r) is plotted as described in 
the text and error bars are suppressed for clarity. 

The main source of error is due to the stochastic noise 
when computing the all-to-all propagators. By imple- 
menting the one-end trick, the four-point function is au- 
tomatically summed over sink and source coordinates and 
thus this method is expected to suppress stochastic noise 
considerably. 

In Fig. [4] we show the pion correlator computed us- 
ing the one-end trick and the direct method as a func- 
tion of the distance from the origin. To avoid having to 
display all lattice points in the graph we replace points 
lying within a cell of size 0.015 fmx0.05 by their aver- 
age. We normalize the correlator by dividing by its value 
at the origin. The errors in Fig. |4] are not shown for 
clarity. As can be seen, we find that the two methods 
yield consistent results for the correlators. This demon- 
strates that excited states are sufficiently suppressed with 
a sink-source separation of 10 time slices. However, at a 
given distance r, the correlator computed using the direct 
method shows more spread than the one computed us- 
ing the one-end trick. That this refiects larger statistical 
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FIG. 5: Comparison between the relative error of the corre- 
lator computed with the one-end trick (blue crosses) and the 
direct method (red circles). 



noise is shown in Fig. [s] where we compare the relative 
errors of the two binned correlators. As can be seen, at 
large distances the maximum relative error exhibited by 
the one-end trick method is around 4% while for the di- 
rect method exceeds 10%. This is a direct consequence 
of the double sum accomplished with the implementa- 
tion of the one-end trick. In addition, when using the 
one-end trick the density-density correlator of a state of 
spin projection = is symmetric under reflections of 
the spatial coordinates i.e. C{f) ~ C{—f) by construc- 
tion whereas in the direct method it is symmetric only 
statistically. For the = ±1 projections of the vec- 
tor meson we instead have C""'=^"'"^(r) = C™^^^^(— r). 
Because of this symmetry we average over the results 
for the niz = +1 and to^ — —1 spin projections and 
hereby denote this correlator by = ±1. The same 
is done for the spin projections = ±3/2 of the A. 
The reduction of the error by more than a factor two 
when using the one-end trick comes at a reduced com- 
putational cost. In the one-end trick the computation of 
the correlator is done using 64 inversions while for the di- 
rect method used in this comparison we carried out 300 
inversions per configuration i.e. we need 4.7 times less 
inversions for twice the accuracy. This, combined with 
the fact that the computation using the one-end trick is 
carried out for a source-sink separation of 14 time slices 
while for the direct method we used a separation of 10 
time slices and given that relative errors grow exponen- 
tially with the sink-source separation, clearly shows the 
superiority of the one-end trick. 

One of the main goals of this calculation is to detect 
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FIG. 6: Comparison between the vector meson — cor- 
relator projected along the three axes computed with the di- 
rect method (upper graph) and with the one-end trick (lower 
graph) using 200 configurations. 



a possible asymmetry in the hadron charge distribution. 
In Fig. [6]we compare the two methods for the case of the 
rriz = spin projection of the vector meson at the lowest 
pion mass available using the same number of configura- 
tions. Only the profile of the correlator along the three 
axes is plotted so that we can detect a possible asymme- 
try. As can be seen, an elongation along the z axis is 
clearly observed only when using the one-end trick. The 
statistical error in the direct method is not small enough 
to draw definite conclusions, since the projections of the 
correlator on the three axes are within error bars. Us- 
ing the one-end trick the fiuctuations are small enough to 
conclude that the vector meson is indeed prolonged along 
the 2 axis. When discussing results on baryon deforma- 
tion one has to keep in mind that statistical fiuctuations 
are larger than for mesons and that we can only apply the 
direct method making reaching conclusions for baryons 
more difficult. 

Having demonstrated the effectiveness of the one-end 
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FIG. 7: Projections of the correlator along the three axes. 
From top to bottom: For the pion and the p - meson us- 
ing the one-end trick with 64 inversions, for the nucleon and 
the A using the direct method with 300 inversion needed per 
configuration. 



B. Results without volume corrections 

In Fig.[7]we show the density-density correlators for the 
pion and the spin zero projection (m^ = 0) of the p-meson 
using the one-end trick and for the nucleon and spin 
= ±1 projection of the A using the direct method. 
All correlators are projected along the three axes to dis- 
play a possible asymmetry. This is done for the smallest 
pion mass available, namely = 0.691(8) GeV. As can 
be seen, a clear elongation of the vector meson along the 
z axis is observed confirming our previous results [16]. 
The asymmetry is clearly smaller than for the lightest 
pion mass shown in Fig. [6] showing that the deformation 
increases as the pion mass decreases. On the other hand, 
the nucleon shows no asymmetry within this method. 
For the A although there is a tendency for results pro- 
jected along the z-axis to lie lower, all projections are 
well within error bars and therefore no asymmetry can 
be claimed. As pointed out when discussing results on 
the p using the direct method, statistical errors can hide 
possible deformation and one may have to improve on 
the errors to detect a small asymmetry. 

Another way to visualize the asymmetry is to construct 
two-dimensional contour plots. Fig. |8] shows a contour 
plot of the rUz = spin /fj-meson state on the x - z plane. 
As can be seen, the contours are elongated along the z- 
axis as compare to a circle of radius equal to the distance 
along the a;-axis for all three pion masses showing a clear 
asymmetry. This leads to the conclusion that the vector 
meson in the spin projection zero state is prolate. On the 
other hand, the rUz = ±1 p-meson state, shown in Fig. [9] 
shows the opposite behavior. Namely the correlator is 
found to be larger along the x-axis, as compared to a 
circle, evidence that in this spin state the p is in fact 
an oblate. This is in agreement with what is derived in 
Ref. [16] where it is shown that if the spin-0 state is a 
prolate the ±1 channels will be oblate with about half the 
amount of deformation. The fact that the p-meson in its 
maximal spin projection state is an oblate is in agreement 
with a recent calculation of a negative electric quadrupole 
form factor evaluated in quenched lattice QCD [51] . 



C. Results after finite volume corrections 

Density-density correlators computed in a finite box 
with periodic boundary conditions are susceptible to fi- 
nite volume effects. Finite volume effects mostly affect 
the tail of the distributions and need to be corrected. To 
perform these corrections we follow the analysis devel- 
oped in Ref. [H]. The density-density correlation func- 
tion computed on a lattice of spatial dimension L can be 
written as an infinite sum over the Brillouin zones 



trick in suppressing stochastic noise, all meson observ- 
ables that we present hereon are computed with the one- 
end trick. 



(20) 



n=Q 
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where C(r) is the density-density correlator computed 
on the periodic lattice and Co{'r) is the "correct" correla- 
tor that one would compute if the lattice were of infinite 
size. Thus the correlation function computed in a finite 
box with periodic boundary conditions is in fact a sum 
of all images arising from the surrounding boxes. Since 
Co(r) is a fast decaying function, approximated by expo- 
nential or Gaussian dependence on the radius, it means 
that the leading contributions to the sum come from the 
nearest neighboring Brillouin zones. A two-dimensional 
sketch drawn in Fig. [TO] demonstrates the images that 
contribute to the correlator. In this figure, the aster- 
isk shows the origin of the fundamental cell (white box) 
while the triangles show the origins of the neighboring 
cells (gray boxes). To first order, the correlator computed 
in the white box is a superposition of the correlator with 
origin the asterisk and the eight correlators with origins 
the filled triangles, in accord with the expression given in 
Eq. ( 20 ) . Thus the correlator that we compute on a peri- 



odic lattice is overestimated. This is particularly severe 
close to the boundaries of the lattice where contributions 
from the images are largest. For example, the correlator 
at the distances indicated by the filled circles in Fig.[TO|is 
approximately twice as large as the "correct" correlator 
since besides the contribution from the fundamental cell, 
a neighboring cell contributes equally as indicated by the 
dashed line. Similarly, the correlator computed at the 
distances indicated by the open circles at the corners of 
the fundamental cell is approximately four times larger 
since there are contributions from three neighboring cells, 
as shown by the dotted line. 

This analysis can be extended to three dimensions. 
The correlator is twice as large at the six distances given 
by ±L/2ni,i = x,y,z where fii is the unit vector in the 
i-direction. Similarly, the correlator is four times as large 
at the twelve distances L/2{hi±nj), i ^ j and eight times 
as large at the eight corners L/2{±.fix ± hy ± hz). 
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FIG. 8: The correlator of the p - meson, — projected on the x - z plane for decreasing pion mass left to right. The dashed 
circles are to guide the eye. 



All results that have been discussed so far are for the 
correlators computed on the lattice with no corrections 
applied for the images. For the analysis of quantities, 
such as the root mean squared radius, that are sensitive 
to the long distance behavior of the distributions it is 
important to take in to account the image contributions 
and define a corrected correlator. To correct for the im- 



ages and extract Cq (f) of Eq. ( 20 1 by knowing only C (f) 



we need to have an Ansatz for the asymptotic behavior of 
Co{r). If the asymptotic behavior is known then we can 
subtract from the lattice data the contribution from the 
images, up to a given order, and extract Co{r). In this 



work, we consider only nearest neighbor contributions to 
the correlator. Thus Eq. (20 1 becomes: 



C{r) 



E Coir- 

|n|G[0,V3] 



nL). 



(21) 



We make an Ansatz for the functional form of Co(r) 
that provides a good description of the data. For instance 
for the pion correlator that is found to be independent of 
the angles, a spherically symmetric Ansatz is used. We 
then perform a least squares fit to the lattice data of the 



sum given on the right hand side of Eq. (21) extracting 
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FIG. 9: The correlator of the p - meson, 
dashed circles are to guide the eye. 



= ±1 projected on the x - z plane for decreasing pion mass left to right. The 
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TABLE II: The parameters obtained from fitting the sum of 
images to the lattice data. 



FIG. 10: Two-dimensional example of image contributions. 
The correlator computed at the filled circles (open circles) 
is approximately two (four) times larger than the "correct" 
correlator. 



the fit parameters of the function that describes Co(r). 
The corrected correlator is then constructed by subtract- 
ing from the lattice data the images determined from the 
fitted function to obtain: 



K 


0.1575 


0.1580 


0.15825 


Mesons 



Ao 0. 

mo 0. 

cr 0. 

Ao 0. 

mo 0. 

Ai 0. 

mi 0. 

cr 1. 

Ao 0. 

mo 0. 

Ai -0. 

mi 0. 

cr 1. 



986(21) 

307(7) 

993(7) 

969(13) 

0173(19) 

00170(31) 

0466(87) 

615(41) 



p, m 



p, m 



,976(10) 

,0194(16) 

,00113(18) 

,0560(91) 

,577(30) 



1.129(33) 
0.405(11) 
0.884(9) 
= 
0.964(21) 
0.0140(26) 
0.0031(16) 
0.077(33) 
1.646(69) 
= ±1 
0.961(16) 
0.0128(16) 
0.00054(34) 
0.025(12) 
1.659(47) 



1.437(78) 
0.579(25) 
0.779(12) 

0.919(31) 

0.0093(26) 

0.00183(46) 

0.0033(12) 

1.76(11) 

0.977(28) 

0.0141(34) 

■0.0012(17) 

0.066(69) 

1.613(87) 



Ao 1 
mo 
cr 1 

Ao 1 
mo 
Ai -0 
mi 
a 1 



.014(39) 

.0673(40) 

.451(20) 

A 

.024(22) 

.0125(11) 

.00029(25) 

.024(13) 

.750(32) 



Baryons 
iV 

1.039(34) 

0.0698(44) 

1.413(22) 

, = ±1 
1.033(19) 
0.0130(12) 

-0.0007(14) 
0.022(25) 
1.708(34) 



1.057(34) 

0.0548(38) 

1.450(24) 

1.023(16) 
0.0087(8) 
-0.00121(49) 
0.077(30) 
1.787(33) 



C'°"{r}^C{r}- Co{f+nL). (22) 

|n|6(0,%/3] 



The Ansatze for Co(r) for each particle are summarized 
below: 
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Corrected 
Uncorrected x 
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r(fm) 

FIG. 11: The pion correlator (top) and the nucleon correlator 
(bottom) as computed on the lattice (crosses) and corrected 
for the images of nearest neighboring lattices (open circles). 
The corrected correlator is divided by a factor of ten for clar- 
ity. Data are binned and error bars are omitted to avoid 
cluttering. 



Q = Ao exp (-TOor" 



exp {—rrii^ir") + Ai exp {-~mir")r^ P2{cos I 



N 



same as for tt, 



= same as for p. 



(23) 



As can be seen, for the pion and the nucleon we 
take spherical functions. For the case of the p we have 
parametrized the correlator in such a way so that an 
asymmetry, as seen in the uncorrected data, is allowed. 
For the A, although no asymmetry can be seen within 
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FIG. 12: The p-meson, = correlator (top) and the A, 
rriz = ±3/2 correlator (bottom). The notation is the same as 
that of Fig. [11] 



our statistical errors we use the same Ansatz as for the 
p to see if the data allow for such a term. 

Since the spatial part of the correlators is even under 
reflection, only L — and L ^ 2 angular momentum 
quantum numbers are allowed. Thus for the p-mcson 
and the A we include an i = 2 component by includ- 
ing the Legendre polynomial P2{cos6). In Table [ll| we 
summarize the fit parameters obtained. The fact that 
for the spin projection — Op state the asymmetric 
term with coefficient Ai is found non-zero and positive 
confirms that the correlator is indeed elongated along the 
2;-axis (prolate) while the same parameter is consistently 
negative for the ruz = ±1 channels pointing to a corre- 
lator larger at the equator (oblate). For the A the Ai 
coefficient comes out negative for all quark masses albeit 
with a large statistical error not allowing any definite 
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conclusions on the A shape. 



In Figs. 11 and 12 we show a comparison between the 
raw lattice data and the lattice data after subtracting 
image contributions for the heaviest pion mass available. 
As can be seen, the correction procedure clearly compen- 



sates for the images, i.e. the spikes at ^/2L/2 and 
are corrected for, leading to a smoother correlator 
that decreases more rapidly at the tails. 



TABLE III: {x^ + y^)/2, {z^) and their difference for each particle at all three pion masses in fm^, left for mesons and right 
for baryons. All errors are jack - knife errors. 



mi (GeV^) (x^ + y')/2 






m'i (GeV^) 








0.477 
0.259 
0.147 


0.1449(6) 
0.1542(7) 
0.1529(7) 


TV 

0.1460(7) 
0.1531(9) 
0.1533(14) 


0.0011(8) 
-0.0010(10) 
0.0005(18) 


0.477 
0.259 
0.147 


0.164(1) 
0.170(1) 
0.181(1) 


N 

0.159(1) 
0.168(2) 
0.182(2) 


-0.006(2) 
-0.002(3) 
0.0008(31) 


0.477 
0.259 
0.147 


P, 

0.174(2) 
0.188(4) 
0.190(5) 


niz — 
0.192(2) 
0.196(6) 
0.207(6) 


0.018(3) 
0.007(7) 
0.016(7) 


0.477 
0.259 
0.147 


A, = ±1 
0.177(1) 0.172(1) 
0.182(1) 0.180(2) 
0.195(2) 0.198(3) 


-0.005(2) 
-0.001(2) 
0.003(4) 


0.477 
0.259 
0.147 


P, 

0.183(1) 
0.199(2) 
0.200(4) 


rriz — ±1 
0.173(2) 
0.186(2) 
0.193(5) 


-0.009(2) 
-0.013(2) 
-0.007(6) 











z fm 




z fm 





0.51 



0.5 




FIG. 13: Three-dimensional contour plot of the p-meson, 
ruz = correlator (red or darker surface) compared to 
a sphere (green or lighter surface). The sphere radius is 
approximately 0.5 fm. The contour shows all r such that 
Cir) = iC(0). 



FIG. 14: Three-dimensional contour plot of the p-meson, 
rriz = ±1 correlator (red or darker surface) compared to 
a sphere (green or lighter surface). The sphere radius is 
approximately 0.5 fm. The contour shows all r such that 
C{r} = |C(0). 



Having corrected the data for the nearest images we 
can now proceed to a quantitati ve a nalysis of the particle 
charge distributions. In Table III we give {x'^ + y^)/2. 



(z^) and their difference for each particle at each pion 
mass available. All errors are jack - knife errors. Here, 
the moments presented are computed using the corrected 
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correlator: 



^.C— (r) ■ 



(24) 



From the difference 



we see once again 



that the = projection of the p is larger along the 
z axis while the = ±1 projections are larger along 
the equator. An additional observation here is that the 
asymmetry of the = ±1 projections is approximately 
half that of the — projection, thus verifying the 
result reached in Ref. [16 . For the case of the A on the 
other hand a spherical distribution cannot be excluded, 
although for the two lightest pion masses we increase 
the statistics by computing the correlators using the first 
and the second half temporal extent of the lattice and 
by using N^—Q noise vectors for the smallest of the two 
values. 

The asymmetry in the p is nicely represented by a 
three-dimensional contour plot. In Figs. 13 and [14] we 
show contour surfaces for the /o-meson in the niz = and 
TOj = ±1 channels respectively, at the intermediate pion 
mass. The correlator is compared to a sphere centered 
at the origin. Once again we see that the = state 
is elongated along the poles while the = ±1 channels 
are flatter. 



VI. SUMMARY AND CONCLUSIONS 

In this work we develop the formalism for the exact 
evaluation of the equal time density-density correlators, 
which in the non-relativistic limit reduce to the hadron 
charge distribution. The pion, p-meson, nucleon and A 
density-density correlators are evaluated using dynami- 
cal Wilson fermions down to a pion mass of 384 MeV. 
The all-to-all propagators needed for the calculation of 
these correlators are computed using stochastic tech- 
niques combined with dilution. Having the all-to-all 
propagators is required so that an explicit projection to 



zero momentum initial and final states is carried out. 
In the meson-sector we implemented the one-end trick, 
which leads to a significant improvement in the accuracy 
with which the density-density correlators are obtained. 
This improved accuracy is needed to conclude with cer- 
tainty that the the p-meson is deformed. The p is found 
to be a prolate when in the spin projection zero state and 
an oblate in the spin projection ±1 state. This result 
corroborates previous studies where the density-density 
correlator of the p was calculated without explicit zero- 
momentum projection and with less accuracy |16| . It is 
also in agreement with a negative quadrupole form factor 
calculated recently on the lattice (34] • For the baryons 
a spherical distribution can not be excluded given the 
present statistical errors despite increase in statistics. 

Finite spatial volume effects affect mainly the long dis- 
tance behavior of the correlators. By adopting an Ansatz 
for the asymptotic dependence of the correlators we cor- 
rect for these finite volume effects by subtracting the first 
image contributions. The functional form determined 
from fits to the corrected data confirm a deformed shape 
for the p meson. For the A, although the fits allow for 
a small deformation, the statistical error is too large to 
exclude a spherical distribution. Further improvements 
in the evaluation of all-to-all propagators such as com- 
bination of stochastic techniques and lower eigenmode 
projection are currently being investigated by a number 
of groups with promising results |35j that have potential 
application in the study of baryon density-density corre- 
lators. 
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